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We study equilibrium sequences of close binary systems composed of identical polytropic stars in 
Newtonian gravity. The solving method is a multi-domain spectral method which we have recently 
developed. An improvement is introduced here for accurate computations of binary systems with 
stiff equation of state (7 > 2). The computations are performed for both cases of synchronized and 
irrotational binary systems with adiabatic indices 7 = 3, 2.5, 2.25, 2 and 1.8. It is found that the 
turning points of total energy along a constant-mass sequence appear only for 7 > 1.8 for synchro- 
■ nized binary systems and 7 > 2.3 for irrotational ones. In the synchronized case, the equilibrium 

' sequences terminate by the contact between the two stars. On the other hand, for irrotational 

, binaries, it is found that the sequences terminate at a mass shedding limit which corresponds to a 

' detached configuration. 
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I. INTRODUCTION 



The research for equilibrium figures of binary systems has a long history. The pioneer work was done by Roche 
', who treated a synchronously rotating, incompressible ellipsoid around a gravitating point source. Following this 
■ work, Darwin constructed a synchronized binary system composed of double incompressible ellipsoids (see [|l|). For 
' the non-synchronized case, Aizenman calculated an equilibrium figure of incompressible ellipsoid with internal motion 
orbiting a point source companion In 1980's, the improvement of computer performance made possible to construct 
equilibrium sequences of synchronized binary systems with compressible equation of state in Newtonian gravity 
Since the beginning of 1990's, the studies on binary systems gain some importance in relation to astrophysical 
I sources of gravitational waves. For example, coalescing binary neutron stars are expected to be one of the most 
I [ promising sources of gravitational radiation that could be detected by the interferometric detectors currently in op- 
^ • eration (TAMA300) or under construction (GEO600, LIGO, and VIRGO). Until now, numerous theoretical studies 
rr' have been done in this research field. Among them there are (i) post-Newtonian analytical treatments (e.g. Q), (ii) 
black hole perturbation |^,^, (iii) Newtonian f^-p^ and post- Newtonian fll|-[r^ (semi-)analytical treatments includ- 



ing hydrodynamical effects of stars, (iv) post-Newtonian hydrodynamical computations |15 1^, (v) fully relativistic 
hydrodynamical treatments, pioneered by the works of Oohara & Nakamura (see e.g. 0), and recently developed by 
' Shibata [p8|-pl| , Oohara & Nakamura |^ , and the Neutron Star Grand Challenge group |23|, p4|. (vi) In parallel of the 
dynamical studies in (v), a quasiequilibrium formulation of the problem has been developed p5|-p8[ and successfully 
implemented [p9|-p3[ . The basic assumption underlying the quasiequilibirum calculations is that the timescale of the 
orbit shrinking is larger than that of the orbital revolution in the pre-coalescing state. Consequently the evolution of 
the binary system can be approximated by a succession of exactly circular orbits, hence the name quasiequilibrium. 

The first quasiequilibrium configurations of binary neutron stars i n g eneral relativity have been obtained four 
years ago by Baumgarte et al. Q,^, followed by Marronetti et al. Ig^l- However these computations considered 
synchronized binaries. As far as coalescing binary neutron stars are concerned, this rotation state does not correspond 
to physical situations, since it has been shown that the gravitational-radiation driven evolution is too rapid for the 
viscous forces to synchronize the spin of each neutron star with the orbit p7| , p8[ | as they do for ordinary stellar binaries. 
Rather, the viscosity is negligible and the fluid velocity circulation (with respect to some inertial frame) is conserved 
in these systems. Provided that the initial spins are not in the millisecond regime, this means that close binary 
configurations are well approximated by zero vorticity (i.e. irrotational) states. Irrotational configurations are more 
complicated to obtain because the fluid velocity does not vanish in the co-orbiting frame (as it does for synchronized 
binaries). 

We have successfully developed a numerical method to tackle this problem and presented the first quasiequilibrium 
configurations of irrotational binary neutron stars elsewhere [ p9| , ^ . The numerical technique relies on a multi-domain 
spectral method within spherical coordinates. Since then, two other groups have obtained relativistic irrotational 
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configurations: (i) Marronetti, Mathews & Wilson [pO|,M by means of single-domain finite difference method within 
Cartesian coordinates and (ii) Uryu & Eriguchi ||4l| , pi| , p2[ by means of multi-domain finite difference method within 
spherical coordinates. 

Recently, we have presented our method in detail with numerous tests |3^ ] (hereafter Paper I). Using this method, 
we have calculated (quasi)equilibrium sequences of binary systems with synchronized and irrotational rotation states. 
In the present article, we will show the results in Newtonian gravity. The results of relativistic calculations will be 
given in a forthcoming article 

The plan of the article is as follows. We start by presenting the equations governing binary systems in Newtonian 
gravity in Sec. ||. Section III is devoted to the brief explanation about the improvement on the cases of stiff equation 



of state such as 7 > 2. In Sec. IV some tests for numerical code are performed. Then, we will show the results of 
evolutionary sequences of both synchronized and irrotational binary systems constructed by double identical stars 
with polytropic equation of state in Sec. 0. Section Vl contains the summary. 



Throughout the present article, we use units of G = c = 1 where G and c denote the gravitational constant and 
speed of light, respectively. 

II. FORMULATION 

Since the method which we use in the present article for getting equilibrium configurations has already been 
explained in Paper I | ^3| , we will only briefiy mention the basic equations and the solving procedure in this section. 

A. Basic equations 

The basic equations governing the problem are as follows (see also Sec. II. F of Paper I): 

1. Equation of state: 

We adopt a simple equation, i.e. a polytropic equation of state; 

p{H) = nn{H)\ (1) 

p being the fluid pressure, n the fluid baryon number density, k a constant, 7 the adiabatic index, and H the 
specific enthalpy. The relation between H and n is given by 



n{H) = 



7 - 1 "IB ^' 

7 K 



1/(7-1) 



(2) 



where me denotes the mean baryon mass (me = 1.66 x 10 kg). 
2. Euler equation: 

dw - 1 - ^ 

— +vVv = Vp~Vu, (3) 

01 niB n 

V being the gravitational potential, and v the velocity field in the inertial frame which is expressed as 

V = u + r2 X r (4) 

in terms of the velocity field in the corotating frame u and the orbital motion r2 x r. In the rigidly rotating 
case, since there is no motion in the corotating frame (u = 0), we obtain 

V = Jl X r. (5) 

On the other hand, in the irrotational case, the velocity field in the inertial frame is potential, i.e. 

V = W, (6) 

where ^ is a scalar function. 
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Under the stationary condition the Euler equation (||) can be integrated as 

1 



+ -^(ri X r)2 = const (7) 



i? + t^+i(V^')^- (r2 xr) •V*^ const (8) 



in the synchronized case and 

in the irrotational case. 

3. Equation of continuity: 

^ + V.(nv)=0. (9) 

This equation is trivially satisfied by stationary rigid rotation. In the irrotational case, we can rewrite Eq. (|^) 
as 

nA* + Vn • V* = (J7 X r) • Vn. (10) 

4. Poisson equation for the gravitational field: 

t^v = AnmBn. (11) 



B. Solving procedure 

(a) First of all, we prepare the equilibrium figures of two spherical stars. In the present computation, we treat 
binary systems composed of equal mass stars with the same equation of state. However, we symmetrize only 
with respect to the orbital plane. Although this treatment elongates the compute time, we do not symmetrize 
the stars because we are planning to study binary systems composed of different mass stars in the series of this 
research. 

(b) The separation between the centers of the two stars is held fixed. Here we define the center as the point of the 
maximum enthalpy (or equivalently maximum density - see Eq. (|^)). 

(c) By setting the central value of the gradient of enthalpy to be zero, we calculate the orbital angular velocity 
(see Sec. IV. D. 2 of Paper I for details). 

(d) For irrotational configurations, the velocity potential "if is obtained by solving Eq. ([l0|). 

(e) Taking into account the gravitational potential from the companion star and the centrifugal force, we calculate 
the new enthalpy field from Eq. (|^) or (||). 

(f) Using the new enthalpy, we search for the location oi H = 0; This defines the stellar surface since = is 
equivalent to p = (see Eqs. (|^) and (^)). Having determined the new stellar surface, we change the position 
of the inner domain boundary to make it fit with the stellar surface. 

(g) By substituting the new enthalpy in the new domain into Eq. (H), the new baryon density is obtained. 



(h) Inserting the new baryon density into the source term of the Poisson equation (|l l|) , we can get the new gravi- 
tational potential. 

(i) Finally, we compare the new enthalpy field with the old one. If the difference is smaller than some fixed threshold, 
the calculation is considered as having converged. If not, return to step (c), and continue. 

During the calculation, we make the baryon mass to converge to a fixed value by varying the central enthalpy. 
Moreover, some relaxation is performed on the gravitational potential, the centrifugal force, and the orbital angular 
velocity, when such a treatment helps the iteration to converge. The explicit procedures are detailed in Paper I. 
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III. IMPROVEMENT ON THE CASES OF STIFF EQUATION OF STATE 



A. Regularization technique 



The numerical method which we use in this series of research [||,|9|j3| is a multi-domain spectral method. In 
general, spectral methods lose much of their accuracy when non-smooth functions are treated because of the so-called 
Gibbs phenomenon. This phenomenon is well known from the most familiar spectral method, namely, the theory of 
Fourier series: the Fourier coefficients (c„) of a function / which is of class but not C^"*"^ decrease as only 
with increasing n. In particular, if the function has some discontinuity, its approximation by a Fourier series does not 
converge towards / at the discontinuity point. The multi-domain spectral method circumvents the Gibbs phenomenon 
p9| . The basic idea is to divide the space into domains chosen so that the physical discontinuities arc located onto 
the boundaries between the domains. 

However, even if the circumvention of the Gibbs phenomenon has been done by introducing mult i- domains, there 
remains the Gibbs phenomenon at the boundaries between the domains when some physical field has infinite derivative 
at the boundary. For example, when we consider a star, its surface coincides with the boundary of a domain. If the 
star has a stiff equation of state, i.e. 7 > 2 in terms of the adiabatic index, the density decreases so rapidly that 
dn/dr diverges at the surface. In order to recover high accuracy, we have developed a method to regularize the density 
profile by extracting its diverging part (see Sec. IV of Ref. Q). 

Here we briefly summarize this method. For a polytrope with adiabatic index 7, the matter density n behaves as 
^1/(7-1) ^ggg gq^ @))- In the steady state configuration, H is Taylor expandable at the neighborhood of the stellar 
surface because the gravitational potential v is Taylor expandable there (cf. Eq. (|^) or (||)). Therefore H vanishes as 
r — R{9, if), where R{9, <p) is the equation of the stellar surface. Consequently n behaves as n ~ [r — R{9, ip)]^/'^'^~^\ 

We introduce a known potential $div such that ridiv ■= A$div has the same pathological behavior as n and such 
that n — TT-div is a regular function (at least more regular than n) . We then numerically solve 

A<I>rcgu = n- ndiv, (12) 
where ^rcgu $ — 'I'div. Consider, for instance, 

<i>div = ne,e,^) (i-e')^"+'\ (13) 

where a — l/(7 ~ 1)- -F is an arbitrary regular function and ^ is a new radial variable such that ^ = 1 at the surface 
of the star |^^. It is easy to see that A^div has a term vanishing at the surface with the same pathological behavior 
as n, i.e., (1 — ^^)". We have indeed 

A$div = AF(1 - e2)("+2) 

-4(a + 2)e(l-^')("+'^aeF 

+ {a + 2)[-6(l - + 4(a + 1)^^(1 - ^TW, (14) 

where A is the Laplacian computed with respect to (^, 6, (p). 

For calculational simplicity, we choose an harmonic function for F(^, 9, ip). Then, we can write $div = X); m '^im^im, 
where 

<fim ■.^e{i-eY"^'^Yr{0,^), (15) 

aim being some numerical coefficients to be determined and Y"^ the standard spherical harmonics. We then obtain 

ndiv = 5^a,™G(0rr(e,^), (16) 

l,m 

with 

Q(e) = (« + 2)[-(4; + 6)(l-a'"+'^e' 

+4(a + l)e('+2)(l-e')"]. (17) 

We now have to determine the values of the coefficients a/^ which give the most regular function rij-egu '■= n — rtdiv 
The criterion which seems to give the best results is the following one. We expand both n and ndiv as truncated series 
of spherical harmonics Yi"^{6, ip) and Chebyshev polynomial Ti(^), 
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I,L,M 

ni^,e,ip)^ J2 naMOYnO,^). (18) 

■i,/,m— 
I,L,M 

ndiv(e,e,^)= E cii^Cum)Yr{e,^), (19) 



where we expand each of the function Ci{£) m a Chcbyshev series as 



ci{£.)^Y.^Hm)- (20) 

1=0 

The values of aim is computed in such a way that the /th coefBcient of the truncated series of rij-egu vanishes: 

aim = — ■ (^1) 

By means of the above procedure, we ehminate in n the pathological term which vanishes as (1 — ^^)" but we 
introduce another pathological term proportional to (1 — ^^)"^^. However, the divergence occurs in a higher order 
derivative of this term so that it has a much weaker effect on the accuracy of the results. The method can be improved 
by taking 

$div = F(t 6, ^)(1 - e')("+''[ai + a2(l - $') 

+a^{l~ef + --- + aK{l-e)''-'] (22) 

instead of Eq. (p^. The coefficients aK are chosen in such a way that the ffist, second, Kth derivatives of riicgu 
vanish at ^ = 1. Let us call K the regularization degree of the procedure. 

The explicit procedure to determine the coefficients ax is as follows: Let us consider the case K — 2 and choose 
an harmonic function for F{S^, 9, ip) for example. In this case, we need two diverging terms, ridivifc (1 < ^ < K), where 
we define 

n^w-.k-.^Y.^'irnChOYnO,^), (23) 

with 

(0 ^{a + k + + 6)(1 - 

+i{a + k)S,''^+^^l-f)'^+''-^]. (24) 

First of all, we expand n, ndiv:i and ndiv:2 by Yi"^{6,(p) and Ti{£_). Then, we can express Cj'{^) as an expansion in a 
Chebyshev series, 

Q'(e) = E^''^^(^)- (25) 

Finally, we solve the simultaneous equations 

nn,n = a\mCli + af„^Clj, (26) 
ni-iim = a}„iClj_i + af„^Cfj_i, (27) 

so that the /th and (/ — l)th coefficients of the truncated series of nrogu vanish. 



B. Illustration and test 



To validate the above regularization technique, we investigate the improvement in the relative error in the virial 
theorem for a spherical static star, which is defined as 
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\W + 3PI 

Virial error = ' ^^j^ (28) 

W being the gravitational potential energy and P the volume integral of the fluid pressure. We give the values of 
the virial error in Fig. |^. In this figure we show the cases of (1) 7 = 3, (2) 7 = 2.5, (3) 7 — 2.25, and (4) 7 = 2, 
and in each small figure except for the case of 7 = 2 for which no Gibbs phenomenon occurs, we have drawn four 
lines: "no regularization" , regularization degree i^T = 1, 2, and 3. It is found from these plots that the regularization 
dramatically improves the accuracy. However, the improvement by the regularization saturates at K — 2. Note here 
that since there is no Gibbs phenomenon in the calculation of a spherical star in the case of 7 = 2, the virial error 
decreases rapidly (exponential decay - evanescent error) with the increase of the number of radial collocation point 
Nr, and reaches the minimum value due to the finite number of digits (15=double precision) used in the numerical 
computation. 

Considering the saturation of the virial error, we will choose the regularization degree K — 2 m the following 
numerical computations for 7 > 2. We note that since the matter density n behaves as i^i/(T~i), the derivatives of 
n of order higher than the second one diverge at the stellar surface even if 7 < 2. Therefore we will also use the 
regularization {K ~ 2) for the cases of 7 < 2. 
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FIG. 1. Relative error in the virial theorem for a spherical star as a function of the number of radial collocation points Nr. 
Solid line with open circle, dotted with open square, dashed with open diamond, and long-dashed with open triangle denote 
the results without regularization, with regularization K = 1, 2, and 3, respectively. 



IV. TESTS OF THE NUMERICAL CODE 



Numerous tests have been already performed in the previous paper in the irrotational case (see Sec. V.B of Paper 
I) . In particular a direct comparison with the numerical results of Uryu & Eriguchi [ p| has been performed and the 
result presented in Table I of Paper I. Here we will focus on some tests in the synchronized case, not presented in 
Paper I. We will also show the relative error on the virial theorem along an evolutionary sequence in both cases, 
extending the test presented in Paper I (Fig. 7) to adiabatic indices different from 2. 



A. Comparison with analytical solutions 

In order to investigate the discrepancy between the results from the numerical code and those from Taniguchi & 
Nakamura's analytic solution |^ in the synchronized case with 7 = 2, we present the relative differences on global 
quantities as functions of the separation in a log- log plot in Fig. |^ (in the same way as Fig. 12 of Paper I for irrotational 
configurations). The relative differences are defined as follows: 
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^num -S^ana /'on\ 

Md2r!Kcp/2' ^ ^ 

^num ^ana /"Qi \ 

O ' ' ' 

"Kcp 

\Spc:num- Spc:ans.\, (32) 

where JIkcp is the Keplerian velocity for point mass particles 

"K.:^(^)"^. (3.) 

In the Taniguchi & Nakamura's analytic solution the global quantities are expanded up to 0[{Ro/d)^]. After 
subtracting these analytic solutions from the numerical one, only the terms of order higher than 0[{Ro / d)^] should 
remain. Indeed, we can see from Fig. |^ that the discrepancies between numerical and analytical solutions for the 
energy E and the relative change in central density Spc are both higher than 0[{Ro/d)^], and those for the angular 
momentum J and the orbital angular velocity fl are both higher than 0[{Ro/dy]. This means that the numerical 
solution and the analytical one agree up to the accuracy of this latter {0[{Ro / d)^]) . 



Relative difference from analytic solution 

Synchronized case (7^2) 

10"' F ' ' ' ' ' ' ' 




Coordinate separation / Radius of a spherical star 



FIG. 2. Relative differences in total energy E, total angular momentum J, orbital angular velocity Q, and relative change in 
central baryon density Spc when comparing the numerical solution with Taniguchi & Nakamura's approximate analytic solution 
along an equilibrium sequence in the synchronized case with j — 2. The horizontal axis denotes logarithmically d/Ro, 
where d is the separation between the two stellar centers, and Ro the stellar radius at infinite separation. The thick solid line 
is a reference one in order to check the inclinations of the results easily. 



B. Virial theorem 



One of the best indicators of the accuracy of numerical solutions for binary systems in Newtonian gravity is the 
relative error in the virial theorem. This error is defined as 

Virial error = , (34) 

where T, W, and P denote respectively the kinetic energy of the binary system, its gravitational potential energy and 
the volume integral of the fluid pressure. By virtue of the virial theorem, the quantity defined by Eq. (§3) should be 
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zero for an exact solution. We show it in Fig. ||. This Figure can be considered as an extension to that presented in 
Paper I (Fig. 7) to adiabatic indices different from 2. It also contains the synchronized case. We can see from Fig. ^ 
that the code is quite accurate as Virial error < 10^^ for dc/Ro > 3 even if 7 = 3, where da denotes the separation 
between two centers of masses of each star. 
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FIG. 3. Relative error in the virial theorem along an evolutionary sequence. The left panel is for synchronized binaries, and 
the right one for irrotational binaries. Solid, dotted, dashed, long-dashed, and dot-dashed lines denote the cases 7 = 3, 2.5, 
2.25, 2, and 1.8, respectively. The horizontal axis denotes dc/Ro, where da is the separation between the centers of masses of 
the two stars, and Ro the stellar radius at infinite separation. 



V. RESULTS 

Using the numerical method explained in the above sections and in Paper I [ p3| , we have constructed equilibrium 
sequences of both synchronized and irrotational binary systems in Newtonian gravity. Here, we consider only the 
case of binary systems composed of stars with equal mass and equation of state. We use 3 domains (one for the fluid 
interior) for each star and the following number of spectral c oeffic ients: Nr x Ng x = 33 x 25 x 24 in each domain. 



A view of the configuration at the energy turning point (Sec. V C below) along a sequence with a 7 = 3 EOS is shown 
in Fig. for synchronized binaries, and in Fig. for irrotational one. We are interested in the turning point of total 
energy (and /or total angular momentum) because it corresponds to the onset of secular instability in the synchronized 
case and that of dynamical instability in the irrotational one at least for the ellipsoidal approximation . 
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FIG. 4. Synchronized binary at the point of minimum energy and angular momentum along a constant-mass sequence 
("last stable orbit"), for an adiabatic index 7 = 3. The grids on the surfaces corresponds to the collocation points in {9, ip) of 
the spectral method. 




FIG. 5. Same as Fig. |j but for an irrotational binary. 



A. Equilibrium sequences 

Our results for equilibrium constant-mass sequences (evolutionary sequences) with adiabatic indices 7 — 
3, 2.5, 2.25, 2 and 1.8 are presented in Tables || and ||. In these tables, d denotes the separation between the 
centers of the two stars. Let us recall that the center of a star is defined as the point of maximum enthalpy (or 
equivalently maximum density). On the other hand, da denotes the orbital separation between centers of masses of 
two stars. Rq, Oi, 02, 03, and ai ^pp are the radius of a spherical star of same mass, the radius parallel to x-axis 
toward the companion star, the radius parallel to y-axis, the radius parallel to z-axis, and the radius parallel to x-axis 
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opposite to the companion star. The {x, y, z) axes are the same as in Fig. 1 of Paper I. pc and pco indicate the central 
density of a star and that of a spherical star of same mass. The normalized quantities 17, J, and E are defined by 

E 

(38) 

where f2, J, and E denote respectively the orbital angular velocity, the total angular momentum, and the total energy, 
and po is the averaged density of a spherical star of the same mass: 



47ri?3' 



Also listed in Tables || and |n| is the ratio 



where {dH/dr)cq,comp [resp. (9-ff/9r)poie] stands for the radial derivative of the enthalpy at the point on the stellar 
surface located in the orbital plane and looking toward the companion star [resp. at the intersection between the 
surface and the axis perpendicular to the orbital plane and going through the stellar center (z axis)]. This quantity is 
useful because the mass shedding limit ("Roche limit") corresponds to % = (cf. Sec. IV. E of Paper I). When x — 0, 
an angular point (cusp) appears at the equator of the star in the direction to the companion. 

An equilibrium sequence terminates either by a contact configuration, which corresponds to d/ai = 2 or by the 
mass shedding limit given by x = 0. Note that these two conditions are not exclusive, as one can have x = at 
contact. By means of our numerical method, it is difficult to get exactly such configurations, so that the final points 
in Tables || and |l| are close to, but not exactly equal to the real end points d/ai = 2 or x = 0. 

For synchronized binaries, the real end point is the contact one In that configuration, the surface of each star 
has a cusp at the contact point. However, our numerical multi-domain method assumes that the boundary between 
each domain is a differentiable surface. Therefore, for computing very close configurations, we stop the adaptation of 
the domain to the surface of the star when x < 0.3 ~ 0.35 (see Sec. IV. E of Paper I). For synchronized binaries, the 
use of adapted domain is not essential, since no boundary condition is set at the stellar surface. The last lines for 
each 7 in Table || are obtained in this way|^. 

On the other hand, in the irrotational case, leaving the adaptation of the domain to the stellar surface, results in 
some numerical error since the solving method for the equation governing the velocity potential [Eq. (|o|)] assumes 
that the domain boundary coincide with the stellar surface (see Appendix B of Paper I). Therefore we stop the 
calculations at some points which are very close to the cusp points but slightly separated. 

The symbol f in Tables || and || indicates the points of minimum total energy (and total angular momentum) along 
the sequence, also called turning points. In the synchronized case, the minimum points exist for 7 > 2, and in the 
irrotational case, they do for 7 > 2.5. In both synchronized and irrotational cases, the minimum points of total energy 
and total angular momentum coincide with each other. It is worth to note here that we cannot exclude the possibility 
of existence of minimum points for 7 < 2 in the synchronized case and for 7 < 2.5 in the irrotational one, because 
we do not calculate up to contact points. However, we suspect that the critical values of the existence of minimum 
points are 7 ^ 1.8 in the synchronized case and 7 ~ 2.3 in the irrotational one. The detailed discussions are given 



later (see Sec. |VQ ) 



We show the total energy, the total angular momentum, the orbital angular velocity, and the relative change in 
central density along a constant mass sequence in Figs. We can see from figures of total energy and total angular 

momentum that the values for synchronized systems are larger than those for irrotational ones. This is because the 
effect of spin of each star in the synchronized case is larger than that in the irrotational one. Such a spin produces 



^Due to the overlapping of the external domains for contact configurations, it was not possible to compute these latter with 
high accuracy. This is why we stopped Table | slightly before d/ai — 2 
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not only the direct effects like spin energy or spin angular momentum but also larger deformations which results in 
larger quadrupole moments. We can see clearly a turning point in energy and angular momentum curves for large 
values of 7 on Figs. ^ - |l This feature will be discussed in Sec. VC . 

It is found from Fig.^ that the central density decreases in both cases of synchronized and irrotational. The 
decrease of the central density in the synchronized case is about one order larger than that in the irrotational one. 
These behaviors are analytically known. For the synchronized case, Chandrasekhar obtained the lowest order change 
in the central density about 70 years ago and Taniguchi & Nakamura have calculated the higher order change 
p3| , and Taniguchi & Nakamura have also shown it for the irrotational case . 

Total energy 
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FIG. 6. Total energy along an evolutionary sequence. The left panel is for synchronized binaries, and the right one for 
irrotational binaries. Solid, dotted, dashed, long-dashed, and dot-dashed lines denote the cases of 7 = 3, 2.5, 2.25, 2, and 1.8, 
respectively. 

Total angular momentum 



Synchronized case Irrotational case 




FIG. 7. Same as Fig. |g but for the total angular momentum. 
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Orbital angular velocity 



Synchronized case 



Irrotational case 



0.4 




I ' ' ' ' ' ' ^ ' ' ' ' ' ^ 

1.2 1.4 1.6 1.8 1.2 1.4 1.6 1.8 

J/(GM'R„)"' J/(GM'R„)"' 



FIG. 8. Orbital angular velocity as a function of total angular momentum along an evolutionary sequence. The lines have 
the same meaning as in Fig. ^. 

Relative change in central density 



Synchronized case Irrotational case 




FIG. 9. Same as Fig. ^ but for the relative change in central baryon density. Note that the two vertical scales are different. 

In Figs. |l^ ~ |l2|, we show isocontours of baryon density in the synchronized case with 7 = 3,2 and 1.8, respectively. 
These figures correspond to the configurations in the last lines (for each adiabatic index 7) in Table |[ Note here that 
the small rough on the stellar surface of 7 = 3 in Fig. HG is an artif act of the graphical software. The solved figure is 



a completely smooth one, thanks to the technique discussed in Sec. [II A 
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Baryon density (z = 0) 



Baryon density {y=0) 





X / R„ 



X / R„ 



FIG. 10. Isocontour of the baryon density of synchronized binaries with 7 = 3 when the separation is d/Ro = 2.941. The 
plots are cross sections of Z — and Y = planes. The thick solid lines denote the stellar surface. The small rough on the 
stellar surface is an artifact of the graphical software. 



Baryon density (z = 0) 



Baryon density (y=0) 






X / Rq 




-2 





X / R 







FIG. 11. Same as Fig. but for 7 = 2 with the separation d/Ro = 2.849. 
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Baryon density (z=0) 



Baryon density (y=0) 





-2 



-2 



FIG. 12. Same as Fig. |l^ but for 7 = 1.8 with the separation d/Ro 



2.828. 



Isocontours of baryon density with 7 = 3, 2 and 1.8 for irrotational binaries are shown in Figs. ^ - |5[ These 
configurations correspond to the semi-final hnes (for each adiabatic index 7) in Table ||. Again note that the small 
rough on the stellar surface of 7 = 3 in Fig. [l^ is an artifact of the graphical software and that the solved figure is a 
completely smooth one, thanks to the technique discussed in Sec. pi A . We do not depict the configurations of the last 
lines of Table || for the following reason. As mentioned above and in Paper I, we make the boundary the inner domain 
fit with the stellar surface. This procedure is essential to accurately solve equilibrium figures in the irrotational case 
(Appendix B of Paper I). However, due to the apparition of the cusp on the stellar surface, we can no longer adapt 
the domain to that surface because all quantities are expressed by summation of a finite number of differentiable 
functions. For very close configurations, just prior to the apparition of the cusp, the surface is highly distorted so that 
there appear unphysical oscillations when using finite series of differentiable functions (Gibbs phenomenon). In the 
last lines for each 7 in Table ||, such oscillations start to be seen, although they do not alter appreciably the global 
quantities. Therefore to avoid any misunderstanding, we choose not to plot the final lines in Table ||, although we 
listed them. 



Baryon density (z = 0) 



3aryon density (y=0) 





-2 

X / R 



-2 









X / R, 
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FIG. 13. Isocontour of the baryon density in the irrotational case with 7 = 3 when the separation is ci/i?o ~ 2.976. The 
plots are cross sections of Z — and Y = planes. The thick solid lines denote the stellar surface. The small rough on the 
stellar surface is an artifact of the graphical software. 



Baryon density (z = 0) 



3aryon density {y=0) 







X / R, 





X / R„ 



FIG. 14. Same as Fig. |l| but for 7 = 2 with the separation d/Ro = 2.917. 



Baryon density (z = 0) 



Baryon density (y=0) 







X / Rq 



-2 





X / R 







FIG. 15. Same as Fig. |l^ but for 7 = 1.8 with the separation d/Ro — 2.932. 

Figures ^ - |l^ show some isocontours of the velocity potential, as well as the velocity field in the co-orbiting frame, 
for the polytropic indices 7 = 3, 2 and 1.8, respectively. In these figures, we show only one of the two star. The 
companion star is located at the position symmetric with respect to the y ^ z plane. The velocity potential ^'o is 
defined as 



vl'o := ^' - Wo • r, 



(41) 
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where Wq is the constant translational velocity field defined as the central value of W := n x r. Note that the vector 
field is tangent to the stellar surface, as it should be. 



psiO (z=0) Velocity w.r.t corotating frame (z = 0) 




-3 -2 -1 -3 -2 -1 

X / R„ X / R„ 



FIG. 16. Contour of velocity potential (left-liand side) and internal velocity field u (right-hand side) with respect to 
the co-orbiting frame in the orbital plane in the irrotational case with 7 = 3 when the separation d/Ro — 2.976. The thick 
solid lines denote the stellar surface. The thin solid and dashed lines in the figure of vefocity potentiai (left-hand side) denote 
positive and negative values, respectively. 



psiO (z=0) Velocity w.r.t corotating frame (z = 0) 




-2 -1 -2 

X / Rq X / R, 

FIG. n. Same as Fig. |f| but for 7 = 2 with d/Ro = 2.9f7. 
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psiO (z=0) 



Velocity w.r.t corotating frame (z = 0) 




B. End points of sequences: contact vs. cusp 



An equilibrium sequence terminates by the contact between the two stars {d/ai = 2) or by a cusp at the onset of 
mass shedding (x = 0). In order to investigate which final fate occurs, it is helpful to display a sequence in the d/ai—x 
plane. This is done in Fig. |l9| where we compare the synchronized sequence with the irrotational one for 7 = 2. It is 
found from this figure that the value of x in the irrotational case is larger than that in the synchronized case for large 
separations. However, when the separation decreases below d/ai 3, x in the irrotational case decreases rapidly 
and becomes smaller than that in the synchronized case for d/ai < 2.5. We magnify the region near the end of the 
sequences in Fig. If we extrapolate the results up to the zero value of x in the figure, we can speculate about the 
final fates of sequences. In the synchronized case, it seems that all the lines will reach x = at d/ai = 2. This means 
that the cusp does not appear before the two stars contact with each other. On the other hand, in the irrotational 
case, it seems that the lines may reach x = before d/ai = 2. The values of d/ai may be d/ai ^ 2.1 for 7 = 3 and 
d/ai ^ 2.25 for 7 = 1.8. This means that the cusp may be created before contact of stars in the irrotational case. 
This behavior in the irrotational case agrees with the results of Uryu . 

It is worth to note here that the point where the cusp may be created is expressed by the orbital separation divided 
by the radius to the companion star (ai). Then, it looks that the cusp appears later along equilibrium sequences for 
higher adiabatic index, i.e., stiffer equation of state. However, if we express the cusp point by the orbital separation 
divided by the radius of a spherical star with the same mass (i?o), the order of cusp appearance will be reversed. 
Although these behaviors are completely different from each other, the origin is the same. The fluid with lower 
adiabatic index is less affected by the tidal force and slightly deformed. It results in the larger d/ai with fixing the 
orbital separation d. 
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7=2 



0.8 



0.6 



0.4 



0.2 



Synchi'onized case 
Irrotational case 



d/a, 

FIG. 19. Equatorial to polar ratio of the radial derivative of enthalpy x as a function of the separation d (normalized by 
the radius ai). The solid and dashed lines denote the cases of synchronized and irrotational fluid states, respectively. 



0.5 



0.4 



0.3 



0.2 



0.1 



Synchronized case 




Irrotational case 




2 2.1 2.2 2.3 2.4 2.5 2 2.1 2.2 2.3 2.4 2.5 
d/a, d/a. 



FIG. 20. Equatorial to polar ratio of the radial derivative of enthalpy x as a function of the separation d (normalized by the 
radius ai). The left (resp. right) panel is for the synchronized (resp. irrotational) binaries. Solid, dotted, dashed, long-dashed, 
and dot-dashed lines denote the cases of 7 = 3, 2.5, 2.25, 2, and 1.8, respectively 



A simple explanation of the difference in x between in the synchronized case and in the irrotational one is as follows. 
Using Eq. ( p] ) and the decomposition of the orbital motion 



we can rewrite Eq. (g) as 



X r = Wo + 



H + v- X rf + ^(V*o - W^)^ = const. 



(42) 



(43) 
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where Wq has been defined in Sec.VA [see Eq. (|4]])] and Ws is the spin part of the orbital motion defined sor star 
No.l as 



W, ■.= n{-yi, XI, 0), 



(44) 



where (xi, yi, zi) are the Cartesian coordinate centered on star 1 (see Fig. 1 in Paper I). In the following explanation, 
we pay particular attention to star 1. Note here that the last term on the left-hand side of E g. is the kinetic 
energy of the velocity field in the corotating frame (see Eq. (P). Comparing Eq. (0) with Eq. (|43[), we see that the 
enthalpy fields for synchronized and irrotational binaries differ precisely by that kinetic energy. 
Then x for each case becomes 



di' 



synch 









) - 






/ cq.comp 





( dv \ ' 
V dzi / pole 

/ oq,comp V 2 / Z OXi 



-{ — ) 

\dzi / pole 



(45) 



(46) 



Of course, due to the difference of deformation of the stars, the gravitational potential v for irrotational binaries is 
different from that for synchronized ones, even if we set two stars at the same orbital separation in the binary systems. 
However, the largest difference between -^^y'^'^^ and x"^"^"* is the last term in the numerator of Eq. (^). We show this 
behavior in the following. First, we divide x iiito two parts; 



synch _ synch 
A Apot 



irrot 



irrot 
Apot 



synch 
' Xfiow ' 

irrot 
A flow 7 



(47) 
(48) 



where 



synch 
Xpot 



synch 
Xfiow 



irrot 
Apot 



irrot 
Afiow 



/ dv 
\dxi 



cq.comp 



,/ d 
\-2 



-(—) 

\dZlJ pole 



0, 



\dxi 



cq,comp 



1 d 



(V^-o 



(—) 

\dzi / pole 



ai 



(for synchronized binaries), 



eq,comp 



'(—) 

V dzi ) pole 



(for irrotational binaries). 



(for irrotational binaries). 



(49) 
(50) 
(51) 

(52) 



Then, we plot the differences (5xpot = Xp™* - Xpot'''' and ^Xflow = Xa™* " Xflow" Fig. y in the 7 = 2 case. 
Here, ^Xpot denotes the difference in x between synchronized and irrotational cases which includes the gravitational 
potential force plus centrifugal force, and (5xfiow denotes one which has relation to the internal flow in the corotating 
frame. 
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0.05 



-0.05 - 



Difference in x 




FIG. 21. Differences in x terms between synchronized and irrotational cases along equilibrium sequences with 7 = 2. 

We can confirm from tliis figure tliat, at small separation, the difference in the gravitational force plus centrifugal 
force is smaller than that in the force related internal flow. Therefore, we will consider only the last term in the 
numerator of Eq. ( ^6| ) below. 

Now, we assume the following form for V^o- 

V^fo = a;i,yi,zi), g{d,xi,yi, zi), h{d,xi,yi, zi)j , (53) 

where /, g and h are some scalar functions of the orbital separation d and the coordinates (xi, yi, zi). Then the last 
term in the numerator of Eq. (46) can be calculated as 



F 



1 d 



(V*o - W, 
df 



2 



[</-»'i;-'--'(ii-o-"lT: 

Accordingly, the surface value at {ri,9i,ifi) = (ai,7r/2,0) <^ (a;i,yi,zi) — (ai,0,0) becomes 



F, 



cq,comp 



li°A^(^ai) + - l) f — (ai) - 1) 

ai dxi \ ai ) \ dx\ ) 



(54) 
(55) 

(56) 



where we have used h{d,xi,yi,0) = because the velocity field is antisymmetric with respect to the x — y plane. 
Note that we simplify the argument list of / and g from (ai,0,0) to (ai). Let us consider the dominant terms in / 
and g which are 



/^A(d)yi, 
g ~ A{d)xi. 



(57) 
(58) 



Here A is a function of the separation d. The forms ( |57| ) and (^8|) can be justified from the studies by Lai, Rasio & 
Shapiro B] or Taniguchi & Nakamura llQ]. Indeed, for ellipsoidal models, A becomes 



A 



ai + ai 



(59) 



and its dependence of d is 0[{d / Ro)~^] because the dominant effect in the deviation from a spherical star is produced 
by the tidal force. From Eqs. (|57)-(|58|) it appears that the first term in the brackets on the right-hand side of Eq. 

:h the second term, so that one is left with 



(p6|) is negligible as compared wit 
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(60) 



Taking into account Eqs. (|58|) and (59), we can have only two types of behaviors for the function g: 



(i) ^ < 1 and |^(«r) < 1, 

ai oxi 

(ii) ^ < 1 and |^(aO > 1. 

ai oxi 

Since g tends to zero for very large orbital separations, the case (i) will occur in the earlier stage of the sequence. In 
this case, fcq.comp becomes negative value so that it makes x"'''°* larger value than ^^y'^'^^ 0. However, if the case (ii) 
actually occur, i^oq.comp can take positive value, and can be smaller than that in the synchronized case. 

In Fig. ^2|, the y-axis component of V^'o, i.e., the function g, and its derivative are shown as a function of the 
coordinate xi (normalized by the surface value ai). This figure corresponds to the last line of the j — 2 case in Table 
||. We can see from this figure that the y-axis component of V^'o remains smaller than that of Ws throughout the 
interior of the star, but the derivative near the stellar surface becomes larger than unity, i.e., the derivative of the 
y-axis component of Wg. This means that the case (ii) really occurs for a very close configuration. It is worth to 
note here that the decrease of x'"°^ by the term F occurs only near the stellar surface. On the contrary, the term F 
increases x"^'°* around the center of the star. 

Next, we show the surface value of y-axis component of V^'o and its derivative along an evolutionary sequence in 
Fig. m. It is found that the y-axis component of V^'o is smaller than that of Ws throughout the sequence even if we 
extrapolate the lines to d/ai ^ 2.1 for 7 = 3 and to d/ai ^ 2.25 for 7 — 1.8. Moreover, the y-axis component of the 
derivative of V^'o becomes larger than that of W^. for every 7 when the orbital separation decreases than d/ai < 2.5. 
This explains why becomes smaller than for d/ai < 2.5 as we have found from Fig. We can also see 



from Fig. 23 that the orbital separation where the derivative of V^'o overcomes that of Wj is larger for smaller 7. 



This fact leads to the earlier appearance of a cusp for smaller 7. 



Y-axis component of grad(*Pg) toward (9, (p)=(7i/2, 0) 

Irrotational case (1^2), d/R„ = 2.851 



^ Q g Y-axis component of grad{f (j)/naj 

' Y-axis component of WyHaj 




Y-axis component of (d(gmd(*i'^})ldx^)/Q. 
Y-axis component of (dW^dXjVn 




x,/a, 



FIG. 22. Y-axis component of V^&o (i.e., the function g) and Ws (upper panel), and their derivatives (lower panel) as a 
function of the coordinate Xi for a close irrotational binary. Solid and dashed lines are the terms concerned with V^'o and Ws, 
respectively. 



^It is worth to note that since the terms —{dv/dxi 
larger for the negative value of F. 



'{du/dzi), and f2^(— d/2 + xi) have negative values, i^'"'™* becomes 
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Y-axis component of gradC^P^) along a sequence 




FIG. 23. Surface value of y-axis component of V^o (upper panel) and its derivative (lower panel) along an evolutionary 
sequence. Solid, dotted, dashed, long-dashed, and dot-dashed lines denote the cases of 7 = 3, 2.5, 2.25, 2, and 1.8, respectively. 
The thick dashed line in the lower panel is the y-axis component of the derivative of Ws . 

Finally, let us discuss the dependence on the orbital separation in the y-axis component of V^o (and also its 
derivative). We present the y-axis component of V^o and its derivative along an evolutionary sequence in log- log 
plot in Fig. |2^. It is found that both y-axis component of V^'o and its derivative behave as 0[{d/ Rq)~^] for the case 
of larger separation jl^. For very close case as cI/Rq < 3, the y-axis component of V^I^o becomes as 0[{d/ Rq)'^"^] 
and its derivative reaches ~ O[(d/i?o)~^^]- 
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Log plot of gradyC^Pp) along a sequence 

Irrotational case (7^2), Values at (r, 6, (p) = (a^, 7t/2, 0) 



grad(>I'„)/na, 
(d(grad (>P„))/dx,)/n 
Line parallel to (d/R^) 
Line parallel to (d/R^) 
Line parallel to (d/R^) 
Line parallel to (d/R^) 



Coordinate separation / Radius of a spherical star 
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FIG. 24. Y-axis component of V^&o and its derivative along an evolutionary sequence. Thick dotted, thick dashed, thick 
long-dashed, and thick dot-dashed lines are reference ones. 
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C. Turning points of total energy 



We show values of x s-nd separations cIq/Rq at which the total energy (and/or total angular momentum) takes its 
minimum along a sequence as a function of the adiabatic index 7 in Fig. ^ We are interested in the turning point 
of total energy (and /or total angular momentum) because it corresponds to the onset of secular instability in the 
synchronized case and that of dynamical instability in the irrotational one at least for the ellipsoidal approxim ation 



Note that the turning points in the total energy and total angular momentum coincide. As discussed in Sec. VB, 
X = denotes the end point of equilibrium sequences for both synchronized and irrotational cases. If we extrapolate 
the results up to the zero value of x, we can obtain the critical value of 7 below which the turning points of the total 
energy does no longer exist, and the value of the corresponding separation. 

We can see from Fig. ^ that x seems to be zero at 7 = 1.7 ~ 1.8 with dc/Ro = 2.7 ^ 2.75 in the synchronized 
case and at 7 = 2.2 ^ 2.3 with dc/Ro = 2.8 ~ 2.85 in the irrotational one. 

Taking into account the appearance the cusp and the existence of the turning point, we can expect that the 
subsequent merger process stars from 

• the turning point for 7 > 1.8 (plunge ?), 

• the contact point for 7 < 1.8 
in the synchronized case, and from 

• the turning point for 7 > 2.3 (plunge ?), 

• the cusp point with mass shedding 7 < 2.3 
in the irrotational one. 



Synchronized case Irrotational case 



- (1) ^^-^^^'^^ 

— 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 


_ (2) 


(3) 


1 1 1 1 1 1 

' (4) 







1.5 1.75 2 2.25 2.5 2.75 2.25 2.5 2.75 3 

Y Y 



FIG. 25. Turning points of total energy as a function of adiabatic index 7. Panels (1) and (3) are for synchronized binaries 
and (2) and (4) are for irrotational ones. 



VI. SUMMARY 



We have studied equilibrium sequences of both synchronized and irrotational binary systems in Newtonian gravity 
with adiabatic indices 7 = 3, 2.5, 2.25, 2 and 1.8. Through the present article, we have understood the qualitative 
differences of physical quantities between synchronized binary systems and irrotational ones. The summary of the 
results is as follows: 

• The two stars contact with each other in the synchronized case as an end point of equilibrium sequence. 
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• Irrotational sequences may terminate instead by a cusp point (mass-shedding) corresponding to a detached 

configuration. 

• The turning points of total energy appear in the cases of 7 > 1.8 in the synchronized case. 

• The turning points of total energy appear in the cases of 7 > 2.3 in the irrotational case. 

• The turning points of total energy and total angular momentum coincide with each other. 
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TABLE I. Orbital angular velocity f2, total angular momentum J, total energy E, axial ratios, relative change in central 
density, equatorial to polar ratio of the radial derivative of enthalpy Xi and virial error along equilibrium sequences in the 
synchronized case, f denotes the minimum points of total energy (and total angular momentum) . 



Synchronized case 



d/Ro 


dc/Ro 


d/ai 


fl 




J 


E 


a2/ai 


aa/ai 


Ol,opp/cil 


(Pc - Pco)/PcO 


X 


Virial 


error 














7 = 3 (n = 


= 0.5) 














8.016 


8.016 


7.970 


7.196 


-2) 


2.043 


-1.172 


0.9940 


0.9904 


0.9993 


-1.065( 


-3) 


0.9838 


5.145( 


-6) 


7.014 


7.014 


6.953 


8.793 


-2) 


1.923 


-1.180 


0.9909 


0.9855 


0.9988 


-1.602 


-3) 


0.9754 


5.054( 


-6) 


6.012 


6.012 


5.927 


1.108 


-1) 


1.798 


-1.191 


0.9852 


0.9767 


0.9977 


-2.582 


-3) 


0.9602 


4.920( 


-6) 


5.010 


5.010 


4.882 


1.458 


-1) 


1.669 


-1.205 


0.9733 


0.9589 


0.9949 


-4.598 


-3) 


0.9288 


4.707( 


-6) 


4.008 


4.007 


3.788 


2.044 


-1) 


1.542 


-1.224 


0.9423 


0.9154 


0.9855 


-9.690( 


-3) 


0.8496 


4.326( 


-6) 


3.507 


3.504 


3.186 


2.512 


-1) 


1.488 


-1.235 


0.9032 


0.8649 


0.9701 


-1.586( 


-2) 


0.7526 


3.994( 


-6) 


t3.116 


3.107 


2.614 


3.043 


-1) 


1.465 


-1.240 


0.8282 


0.7781 


0.9306 


-2.672 


-2) 


0.5652 


3.679( 


-6) 


3.006 


2.991 


2.379 


3.250 


-1) 


1.470 


-1.239 


0.7774 


0.7245 


0.8955 


-3.299 


-2) 


0.4277 


4.940( 


-6) 


2.941 


2.913 


2.050 


3.415 


-1) 


1.485 


-1.230 


0.6807 


0.6299 


0.8050 


-3.971( 


-2) 


0.08238 


2.852( 


-3) 














7 


= 2.5 (n 


= 2/3) 














8.197 


8.197 


8.155 


6.959 


-2) 


2.061 


-1.137 


0.9949 


0.9917 


0.9994 


-1.252( 


-3) 


0.9841 


7.434( 


-7) 


6.968 


6.968 


6.908 


8.881 


-2) 


1.914 


-1.147 


0.9915 


0.9864 


0.9988 


-2.057 


-3) 


0.9738 


7.260( 


-7) 


5.738 


5.738 


5.648 


1.189 


-1) 


1.758 


-1.161 


0.9842 


0.9751 


0.9973 


-3.751 


-3) 


0.9518 


7.012( 


-7) 


4.918 


4.918 


4.791 


1.499 


-1) 


1.650 


-1.173 


0.9738 


0.9597 


0.9947 


-6.111 


-3) 


0.9211 


6.747( 


-7) 


4.099 


4.098 


3.901 


1.974 


-1) 


1.544 


-1.189 


0.9510 


0.9274 


0.9878 


-1.115 


-2) 


0.8553 


6.245( 


-7) 


3.689 


3.688 


3.428 


2.319 


-1) 


1.494 


-1.198 


0.9276 


0.8961 


0.9791 


-1.615 


-2) 


0.7895 


5.889( 


-7) 


t3.037 


3.029 


2.511 


3.161 


-1) 


1.443 


-1.210 


0.8201 


0.7701 


0.9212 


-3.626 


-2) 


0.4906 


8.134( 


-7) 


2.951 


2.938 


2.310 


3.329 


-1) 


1.446 


-1.209 


0.7742 


0.7224 


0.8866 


-4.281 


-2) 


0.3524 


3.171( 


-6) 


2.901 


2.883 


2.081 


3.449 


-1) 


1.455 


-1.207 


0.7073 


0.6568 


0.8237 


-4.837( 


-2) 


0.1154 


2.539( 


-4) 














T 


= 2.25 (n 


= 0.8) 














8.164 


8.164 


8.123 


7.001 


-2) 


2.055 


-1.108 


0.9951 


0.9921 


0.9994 


-1.459 


-3) 


0.9835 


1.180( 


-7) 


6.804 


6.804 


6.743 


9.203 


-2) 


1.891 


-1.119 


0.9913 


0.9861 


0.9987 


-2.545 


-3) 


0.9709 


1.143( 


-7) 


5.783 


5.783 


5.697 


1.175 


-1) 


1.760 


-1.131 


0.9854 


0.9771 


0.9975 


-4.205 


-3) 


0.9516 


1.118( 


-7) 


4.763 


4.762 


4.630 


1.573 


-1) 


1.625 


-1.147 


0.9725 


0.9579 


0.9942 


-7.767( 


-3) 


0.9098 


1.055( 


-7) 


4.082 


4.082 


3.890 


1.985 


-1) 


1.535 


-1.161 


0.9532 


0.9306 


0.9881 


-1.290 


-2) 


0.8494 


9.882( 


-8) 


3.402 


3.400 


3.082 


2.625 


-1) 


1.453 


-1.177 


0.9056 


0.8686 


0.9684 


-2.475 


-2) 


0.7058 


8.189( 


-8) 


3.062 


3.056 


2.585 


3.104 


-1) 


1.427 


-1.183 


0.8422 


0.7949 


0.9327 


-3.848 


-2) 


0.5178 


2.051( 


-7) 


t2.980 


2.972 


2.429 


3.249 


-1) 


1.425 


-1.184 


0.8119 


0.7623 


0.9117 


-4.398( 


-2) 


0.4256 


8.013( 


-7) 


2.878 


2.863 


2.104 


3.463 


-1) 


1.430 


-1.182 


0.7262 


0.6762 


0.8356 


-5.426( 


-2) 


0.1355 


4.468( 


-5) 


8.021 


8.021 


7.980 


7.189 


-2) 


2.035 


-1.061 


7 = 2 (n 
0.9952 


= 1) 
0.9923 


0.9994 


-1.820 


-3) 


0.9819 


3.154( 


-13) 


6.806 


6.806 


6.748 


9.199 


-2) 


1.887 


-1.072 


0.9920 


0.9872 


0.9988 


-3.003 


-3) 


0.9698 


3.234( 


-12) 


5.834 


5.834 


5.753 


1.159 


-1) 


1.762 


-1.083 


0.9869 


0.9794 


0.9977 


-4.826 


-3) 


0.9511 


8.484( 


-13) 


4.861 


4.861 


4.740 


1.525 


-1) 


1.631 


-1.098 


0.9763 


0.9635 


0.9950 


-8.548 


-3) 


0.9126 


1.530( 


-12) 


3.889 


3.889 


3.680 


2.135 


-1) 


1.499 


-1.119 


0.9486 


0.9245 


0.9857 


-1.775 


-2) 


0.8160 


2.063( 


-11) 


3.403 


3.402 


3.098 


2.618 


-1) 


1.439 


-1.131 


0.9134 


0.8788 


0.9704 


-2.848 


-2) 


0.6990 


2.283( 


-9) 


3.014 


3.010 


2.532 


3.169 


-1) 


1.403 


-1.140 


0.8431 


0.7970 


0.9291 


-4.646( 


-2) 


0.4725 


2.259( 


-7) 


t2.892 


2.885 


2.268 


3.394 


-1) 


1.399 


-1.141 




0.7373 


0.8845 


-5.676 


-2) 


0.2862 


4.158( 


-6) 


2.849 


2.839 


2.092 


3.487 


-1) 


1.400 


-1.141 


0.7361 


0.6878 


0.8364 


-0. i / D 


-Z 1 


0.1116 


1.315( 


-4) 














7 


= 1.8 (n = 


= 1.25) 














8.097 


8.097 


8.058 


7.088 


-2) 


2.041 


-0.9942 


0.9957 


0.9931 


0.9994 


-2.087 


-3) 


0.9817 


1.829( 


-9) 


6.701 


6.701 


6.643 


9.416 


-2) 


1.869 


-1.006 


0.9923 


0.9877 


0.9988 


-3.712( 


-3) 


0.9672 


1.790( 


-9) 


5.584 


5.584 


5.498 


1.238 


-1) 


1.722 


-1.020 


0.9862 


0.9784 


0.9974 


-6.502 


-3) 


0.9419 


1.790( 


-9) 


4.746 


4.746 


4.622 


1.580 


-1) 


1.606 


-1.034 


0.9764 


0.9639 


0.9947 


-1.081 


-2) 


0.9023 


1.712( 


-9) 


3.909 


3.908 


3.710 


2.118 


-1) 


1.489 


-1.053 


0.9537 


0.9317 


0.9867 


-2.029 


-2) 


0.8135 


1.586( 


-9) 


3.350 


3.349 


3.043 


2.677 


-1) 


1.415 


-1.068 


0.9152 


0.8818 


0.9692 


-3.462( 


-2) 


0.6719 


4.648( 


-9) 


3.071 


3.069 


2.651 


3.062 


-1) 


1.384 


-1.076 


0.8718 


0.8304 


0.9439 


-4.821( 


-2) 


0.5203 


6.411( 


-8) 


2.932 


2.928 


2.405 


3.296 


-1) 


1.372 


-1.079 


0.8296 


0.7840 


0.9135 


-5.863 


-2) 


0.3764 


9.670( 


-7) 


2.828 


2.822 


2.104 


3.496 


-1) 


1.367 


-1.080 


0.7532 


0.7067 


0.8440 


-6.939 


-2) 


0.1175 


1.744( 


-5) 



26 



TABLE II. Orbital angular velocity f2, total angular momentum J, total energy E, axial ratios, relative change in central 
density, equatorial to polar ratio of the radial derivative of enthalpy Xj and virial error along equilibrium sequences in the 
irrotational case, f denotes the minimum points of total energy (and total angular momentum). 



Irrotational case 



d/Ro 


dc/Ro 


d/ai 







J 


E 






Ol,opp/ffll 


(Pc - Pco)/pcO 


X 


Virial 


error 


8.016 


8.016 


7.983 


7.196 


-2) 


2.002 


7 

-1.173 


= 3 (n = 
0.9940 


0.5) 
0.9940 


0.9993 


-7.249 


-6) 


0.9897 


5.194( 


-6) 


7.014 


7.014 


6.970 


8.793 


-2) 


1.873 


-1.182 


0.9909 


0.9910 


0.9988 


-1.638 


-5) 


0.9843 


5.123( 


-6) 


6.012 


6.012 


5.950 


1.108 


-1) 


1.734 


-1.194 


0.9851 


0.9853 


0.9977 


-4.243 


-5) 


0.9742 


5.029( 


-6) 


5.010 


5.010 


4.914 


1.458 


-1) 


1.584 


-1.211 


0.9729 


0.9736 


0.9948 


-1.337 


-4) 


0.9525 


4.891( 


-6) 


4.008 


4.007 


3.838 


2.042 


-1) 


1.420 


-1.235 


0.9408 


0.9434 


0.9853 


-5.862 


-4) 


0.8939 


4.642( 


-6) 


3.507 


3.505 


3.248 


2.506 


-1) 


1.334 


-1.252 


0.8994 


0.9055 


0.9695 


-1.559 


-3) 


0.8138 


4.463( 


-6) 


3.006 


2.993 


2.442 


3.230 


-1) 


1.260 


-1.270 


0.7622 


0.7795 


0.8892 


-6.832 


-3) 


0.4775 


2.016( 


-7) 


t2.976 


2.960 


2.346 


3.296 


-1) 


1.259 


-1.270 


0.7362 


0.7547 


0.8686 


-7.981 


-3) 


0.3973 


9.366( 


-6) 


2.956 


2.936 


2.257 


3.346 


-1) 


1.260 


-1.270 


0.7099 


0.7292 


0.8458 


-9.095 


-3) 


0.3097 


3.220( 


-5) 


8.197 


8.197 


8.168 


6.959 


-2) 


2.025 


7 = 
-1.138 


2.5 (n = 
0.9948 


= 2/3) 
0.9949 


0.9994 


-7.110 


-6) 


0.9900 


7.487( 


-7) 


6.968 


6.968 


6.926 


8.880 


-2) 


1.867 


-1.149 


0.9914 


0.9915 


0.9988 


-1.919 


-5) 


0.9833 


7.382( 


-7) 


5.738 


5.738 


5.674 


1.189 


-1) 


1.694 


-1.164 


0.9841 


0.9843 


0.9973 


-6.355 


-5) 


0.9686 


7.194( 


-7) 


4.918 


4.918 


4.826 


1.498 


-1) 


1.570 


-1.178 


0.9736 


0.9742 


0.9947 


-1.677 


-4) 


0.9475 


7.037( 


-7) 


4.099 


4.098 


3.952 


1.973 


-1) 


1.435 


-1.199 


0.9502 


0.9521 


0.9878 


-5.524 


-4) 


0.8996 


6.760( 


-7) 


3.689 


3.688 


3.490 


2.316 


-1) 


1.364 


-1.212 


0.9260 


0.9297 


0.9791 


-1.149 


-3) 


0.8481 


6.403( 


-7) 


3.279 


3.276 


2.978 


2.778 


-1) 


1.293 


-1.227 


0.8777 


0.8853 


0.9574 


-2.832 


-3) 


0.7372 


6.161( 


-7) 


2.951 


2.941 


2.401 


3.306 


-1) 


1.244 


-1.240 


0.7680 


0.7833 


0.8862 


-7.947 


-3) 


0.4184 


9.912( 


-6) 


t2.902 


2.887 


2.219 


3.415 


-1) 


1.240 


-1.241 


0.7162 


0.7334 


0.8415 


-1.023 


-2) 


0.2305 


7.755( 


-5) 














7 = 


2.25 (n 


= 0.8) 














8.164 


8.164 


8.137 


7.000 


-2) 


2.021 


-1.109 


0.9951 


0.9951 


0.9994 


-7.639 


-6) 


0.9896 


1.179( 


-7) 


6.804 


6.804 


6.762 


9.203 


-2) 


1.845 


-1.121 


0.9913 


0.9914 


0.9987 


-2.327 


-5) 


0.9814 


1.164( 


-7) 


5.783 


5.783 


5.724 


1.175 


-1) 


1.701 


-1.134 


0.9854 


0.9856 


0.9975 


-6.332 


-5) 


0.9686 


1.142( 


-7) 


4.763 


4.762 


4.669 


1.573 


-1) 


1.545 


-1.152 


0.9723 


0.9730 


0.9942 


-2.145 


-4) 


0.9399 


1.108( 


-7) 


4.082 


4.082 


3.943 


1.984 


-1) 


1.432 


-1.170 


0.9528 


0.9544 


0.9882 


-5.859 


-4) 


0.8960 


1.065( 


-7) 


3.402 


3.400 


3.159 


2.620 


-1) 


1.312 


-1.193 


0.9043 


0.9094 


0.9687 


-2.119 


-3) 


0.7810 


9.513( 


-8) 


3.062 


3.057 


2.681 


3.091 


-1) 


1.254 


-1.208 


0.8401 


0.8501 


0.9333 


-5.056 


-5) 


0.6038 


2.917( 


-7) 


2.926 


2.918 


2.406 


3.334 


-1) 


1.233 


-1.213 


0.7806 


0.7942 


0.8903 


-8.038 


-3) 


0.4036 


1.281( 


-5) 


2.875 


2.864 


2.235 


3.442 


-1) 


1.228 


-1.215 


0.7333 


0.7487 


0.8494 


-1.011 


-2) 


0.2226 


8.615( 


-5) 


8.021 


8.021 


7.995 


7.189 


-2) 


2.003 


7 

-1.062 


= 2 {n-- 
0.9952 


= 1) 
0.9952 


0.9994 


-8.767 


-6) 


0.9886 


3.324( 


-12) 


6.806 


6.806 


6.768 


9.198 


-2) 


1.845 


-1.073 


0.9920 


0.9920 


0.9988 


-2.382 


-5) 


0.9808 


2.828( 


-12) 


5.834 


5.834 


5.780 


1.159 


-1) 


1.708 


-1.086 


0.9869 


0.9871 


0.9977 


-6.136 


-5) 


0.9684 


2.685( 


-13) 


4.861 


4.861 


4.779 


1.525 


-1) 


1.560 


-1.103 


0.9763 


0.9767 


0.9950 


-1.912 


-4) 


0.9422 


2.996( 


-13) 


3.889 


3.889 


3.743 


2.134 


-1) 


1.397 


-1.128 


0.9486 


0.9504 


0.9859 


-8.118 


-4) 


0.8724 


1.558( 


-11) 


3.403 


3.402 


3.182 


2.614 


-1) 


1.311 


-1.146 


0.9138 


0.9178 


0.9712 


-2.058 


-3) 


0.7799 


1.467( 


-9) 


2.965 


2.962 


2.561 


3.240 


-1) 


1.232 


-1.166 


0.8295 


0.8390 


0.9536 


-6.198 


-3) 


0.5118 


2.299( 


-6) 


2.917 


2.912 


2.462 


3.328 


-1) 


1.224 


-1.168 


0.8081 


0.8187 


0.9042 


-7.255 


-3) 


0.4298 


8.920( 


-6) 


2.851 


2.844 


2.278 


3.457 


-1) 


1.214 


-1.171 


0.7602 


0.7730 


0.8628 


-9.293 


-3) 


0.2325 


8.369( 


-5) 


8.097 


8.097 


8.073 


7.088 


-2) 


2.012 


7 = 
-0.9951 


1.8 (n = 
0.9957 


= 1.25) 
0.9958 


0.9994 


-8.152 


-6) 


0.9885 


1.846( 


-9) 


6.701 


6.701 


6.665 


9.416 


-2) 


1.831 


-1.008 


0.9923 


0.9924 


0.9988 


-2.573 


-5) 


0.9791 


1.823( 


-9) 


5.584 


5.584 


5.529 


1.238 


-1) 


1.671 


-1.023 


0.9862 


0.9864 


0.9974 


-7.867 


-5) 


0.9624 


1.868( 


-9) 


4.746 


4.746 


4.666 


1.580 


-1) 


1.541 


-1.039 


0.9766 


0.9770 


0.9948 


-2.162 


-4) 


0.9355 


1.805( 


-9) 


3.909 


3.909 


3.777 


2.116 


-1) 


1.400 


-1.061 


0.9543 


0.9556 


0.9871 


-7.505 


-4) 


0.8722 


1.747( 


-9) 


3.350 


3.350 


3.140 


2.673 


-1) 


1.299 


-1.082 


0.9176 


0.9209 


0.9709 


-2.141 


-3) 


0.7629 


3.631( 


-9) 


3.071 


3.070 


2.775 


3.054 


-1) 


1.247 


-1.095 


0.8777 


0.8834 


0.9482 


-4.076 


-3) 


0.6312 


1.613( 


-7) 


2.932 


2.929 


2.550 


3.284 


-1) 


1.221 


-1.102 


0.8399 


0.8475 


0.9213 


-5.944 


-3) 


0.4843 


3.355( 


-6) 


2.837 


2.833 


2.309 


3.461 


-1) 


1.204 


-1.107 


0.7803 


0.7907 


0.8676 


-8.223 


-3) 


0.2118 


1.018( 


-4) 
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